Esta guía explica cómo realizar algunas tareas comunes de análisis de datos de recuento de secuenciación microbiana utilizando tidytacos. Ilustraremos utilizando un conjunto de datos con muestras de microbioma humano del tracto vaginal, tomadas de este artículo por Lebeer et al.

Tidytacos: la filosofía del paquete

Un objeto tidytacos es simplemente una lista de tres tablas:

El paquete se llama tidytacos porque cada una de las tablas está ordenada (en inglés tidy): cada fila representa una observación y cada columna una variable (puedes encontrar más información sobre la ordenación de datos en este enlace).

enlace a dictionario

Configuración

En caso de que aún no hayas instalado tidytacos, puedes instalarlo usando devtools:

install.packages("devtools")
devtools::install_github("LebeerLab/tidytacos")

Para esta guía, sólo necesitamos cargar tres paquetes: tidytacos (por supuesto), el conjunto de paquetes tidyverse, y VTutorials.

library(tidyverse)
library(tidytacos)
library(VTutorials)

Exploremos el dataset vdata

Nuestro conjunto de datos de ejemplo está disponible en el paquete VTutorials y no es necesario importarlo ni convertirlo. Nuestro set de datos de ejemplo Se llama “vdata”.

“vdata” es un objeto de tipo tidytacos y contiene 3 tablas: “counts”, “samples”, y “taxa”.

summary(vdata)
##         Length Class  Mode
## counts   3     tbl_df list
## samples 28     tbl_df list
## taxa    12     tbl_df list

Comenzamos inspeccionando la tabla de muestras:

glimpse(vdata$samples)
## Rows: 200
## Columns: 28
## $ sample                               <chr> "SAMEA13411869", "SAMEA13410951",…
## $ sample_id                            <chr> "s16", "s100", "s114", "s170", "s…
## $ Technical.run_20201009_isala_cross_1 <lgl> FALSE, TRUE, FALSE, FALSE, TRUE, …
## $ Technical.run_20201016_isala_cross_2 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201023_isala_cross_3 <lgl> TRUE, FALSE, TRUE, FALSE, FALSE, …
## $ Technical.run_20201030_isala_cross_4 <lgl> FALSE, FALSE, FALSE, TRUE, FALSE,…
## $ Technical.run_20201106_isala_cross_5 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201120_isala_cross_6 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201204_isala_cross_8 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20201211_isala_cross_9 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.run_20210119_isala_cross_7 <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ Technical.libsize                    <dbl> 37957, 22782, 32150, 7707, 24123,…
## $ General.Age                          <dbl> 24, 60, 24, 33, 29, 24, 23, 42, 3…
## $ Health.BMI                           <dbl> 30.11938, 19.00391, 25.03414, 23.…
## $ Health.Antibiotic.3months            <lgl> FALSE, TRUE, TRUE, TRUE, FALSE, F…
## $ Sexual.Intercourse.24hours           <lgl> FALSE, FALSE, FALSE, FALSE, FALSE…
## $ plate                                <chr> "PlateBA", NA, "PlateBB", "PlateB…
## $ volume                               <dbl> 2.2025248, NA, 1.5580467, 1.83601…
## $ dna_conc                             <dbl> 11.169, NA, 15.789, 20.648, 24.65…
## $ read_conc                            <dbl> 17233.4038, NA, 20634.8110, 4197.…
## $ read_conc_norm                       <dbl> 0.9693078, NA, 1.1606229, 0.57860…
## $ lib_size                             <dbl> 37957, NA, 32150, 7707, 24123, 29…
## $ `Sample month`                       <dbl> 7, NA, 7, 7, 7, 7, 7, 7, 7, 7, 7,…
## $ `Transit time`                       <dbl> 1, NA, 1, 1, 1, 2, 1, 1, 6, 2, 1,…
## $ valencia_subcst                      <chr> "I-A", "IV-C0", "III-A", "IV-B", …
## $ valencia_cst                         <chr> "I", "IV-C", "III", "IV-B", "IV-B…
## $ max_genus1                           <chr> "L. crispatus", "Prevotella", "L.…
## $ max_genus2                           <chr> "Gardnerella", "Other", "Prevotel…

Continuamos con la tabla “taxa”:

glimpse(vdata$taxa)
## Rows: 942
## Columns: 12
## $ taxon             <chr> "t1005", "t1019", "t1023", "t1025", "t1030", "t1031"…
## $ taxon_id          <chr> "t1005", "t1019", "t1023", "t1025", "t1030", "t1031"…
## $ kingdom           <chr> "Bacteria", "Bacteria", "Bacteria", "Bacteria", "Bac…
## $ phylum            <chr> "Bacteroidetes", "Bacteroidetes", "Bacteroidetes", "…
## $ class             <chr> "Bacteroidia", "Bacteroidia", "Bacteroidia", "Bacter…
## $ order             <chr> "Bacteroidales", "Bacteroidales", "Bacteroidales", "…
## $ family            <chr> "EU845084_f", "Porphyromonadaceae", "EU845084_f", "E…
## $ genus             <chr> NA, "Porphyromonas", "DQ003626_g", NA, "Parabacteroi…
## $ species           <chr> NA, NA, NA, NA, "johnsonii", "merdae", NA, NA, "putr…
## $ sequence          <chr> "TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAG…
## $ genus_oldtaxonomy <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ genus_nosubgroups <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …

Ahora veamos la tabla “counts”:

glimpse(vdata$counts)
## Rows: 6,714
## Columns: 3
## $ count     <dbl> 15, 80, 8, 3, 7, 6, 20, 72, 49, 65, 4, 45, 5, 2, 18, 140, 2,…
## $ sample_id <chr> "s470", "s879", "s894", "s1145", "s1218", "s1388", "s1631", …
## $ taxon_id  <chr> "t1005", "t1005", "t1005", "t1005", "t1005", "t1005", "t1005…

Luego echamos un vistazo rápido al número total de muestras (n_samples), ASV (n_taxa), y lecturas (n_reads) en el objeto tidytacos:

tacosum(vdata)
## n_samples    n_taxa   n_reads 
##       200       942   5022096

Hagamos un diagrama de barras apilado de un subconjunto de muestras

Podemos crear muy fácilmente un gráfico para explorar un subconjunto de nuestras muestras (por ejemplo, solo muestras de participantes que han tomado antibióticos en los últimos 3 meses y son menores de 40 años) de la siguiente manera:

vdata %>%
  filter_samples(Health.Antibiotic.3months == TRUE, General.Age <= 40) %>%
  tacoplot_stack()

La función filter_samples hace lo que dice: filtrar muestras. También eliminará los taxones de la tabla de taxones que tengan cero lecturas totales en las muestras restantes. La función tacoplot_stack devuelve una buena visualización de diagramas de barras apiladas de los taxones más abundantes en nuestras muestras.

¿Se te ocurren otras ideas de como filtrar y rápidamente visualizar vdata?

Subconjunto de datos

Nuestra siguiente pregunta para este conjunto de datos es hasta qué punto la actividad sexual dentro de las 24 horas antes de tomada las muestras afecta la composición microbiana de las participantes en edad reproductiva (asumamos <=40 años). Para hacernos una idea primero lo podemos visualizar:

vdata_less40 <- vdata %>% filter_samples(General.Age <= 40)
tacoplot_stack(vdata_less40)+
  geom_point(aes(y=-0.02,color=Sexual.Intercourse.24hours)) +
  geom_point(aes(y=-0.05,color=valencia_cst))

## Diversidad Alpha Para explorar la diversidad alfa, creemos una versión enrarecida (rarefied en inglés) del conjunto de datos:

vdata_rar <- vdata %>%
  add_total_count() %>%
  filter_samples(total_count >= 2000) %>%
  rarefy(2000) %>%
  add_alpha()

La función add_total_count agregará números totales de lectura de muestra a la tabla de muestra.

La función rarefy submuestreará aleatoriamente todas las muestras n veces. Solo funciona si el recuento de lecturas de cada muestra es igual o superior a n. 

Para determinar la riqueza de ASV, optamos por enrarecer primero, pero esto puede depender de sus datos.

La función add_alpha se puede utilizar para agregar varias métricas de diversidad alfa a la tabla de muestra.

Ahora recalculemos el conteo total de lecturas por cada muestra. Debería ser de 2000 lecturas para todas las muestras.

vdata_rar %>%
  mutate_samples(old_total_count = total_count) %>%
  select_samples(-total_count)%>%
  add_total_count()
## $counts
## # A tibble: 5,055 × 3
##    count sample_id taxon_id
##    <int> <chr>     <chr>   
##  1     2 s470      t1005   
##  2    12 s879      t1005   
##  3     4 s1631     t1005   
##  4     7 s1817     t1005   
##  5     2 s1902     t1005   
##  6     4 s1942     t1005   
##  7     4 s2790     t1005   
##  8     1 s100      t1019   
##  9     1 s1902     t1023   
## 10    10 s1331     t1025   
## # ℹ 5,045 more rows
## 
## $samples
## # A tibble: 200 × 32
##    sample        sample_id Technical.run_20201009_isala…¹ Technical.run_202010…²
##    <chr>         <chr>     <lgl>                          <lgl>                 
##  1 SAMEA13411869 s16       FALSE                          FALSE                 
##  2 SAMEA13410951 s100      TRUE                           FALSE                 
##  3 SAMEA13411960 s114      FALSE                          FALSE                 
##  4 SAMEA13412020 s170      FALSE                          FALSE                 
##  5 SAMEA13410958 s177      TRUE                           FALSE                 
##  6 SAMEA13412032 s183      FALSE                          FALSE                 
##  7 SAMEA13412067 s222      FALSE                          FALSE                 
##  8 SAMEA13412068 s223      FALSE                          FALSE                 
##  9 SAMEA13412130 s288      FALSE                          FALSE                 
## 10 SAMEA13412143 s302      FALSE                          FALSE                 
## # ℹ 190 more rows
## # ℹ abbreviated names: ¹​Technical.run_20201009_isala_cross_1,
## #   ²​Technical.run_20201016_isala_cross_2
## # ℹ 28 more variables: Technical.run_20201023_isala_cross_3 <lgl>,
## #   Technical.run_20201030_isala_cross_4 <lgl>,
## #   Technical.run_20201106_isala_cross_5 <lgl>,
## #   Technical.run_20201120_isala_cross_6 <lgl>, …
## 
## $taxa
## # A tibble: 688 × 12
##    taxon taxon_id kingdom  phylum      class order family genus species sequence
##    <chr> <chr>    <chr>    <chr>       <chr> <chr> <chr>  <chr> <chr>   <chr>   
##  1 t1005 t1005    Bacteria Bacteroide… Bact… Bact… EU845… <NA>  <NA>    TACGGAG…
##  2 t1019 t1019    Bacteria Bacteroide… Bact… Bact… Porph… Porp… <NA>    TACGGAG…
##  3 t1023 t1023    Bacteria Bacteroide… Bact… Bact… EU845… DQ00… <NA>    TACGGAG…
##  4 t1025 t1025    Bacteria Bacteroide… Bact… Bact… EU845… <NA>  <NA>    TACGGAG…
##  5 t1030 t1030    Bacteria Bacteroide… Bact… Bact… Porph… Para… johnso… TACGGAG…
##  6 t1040 t1040    Bacteria Bacteroide… Bact… Bact… Riken… Alis… <NA>    TACGGAG…
##  7 t1045 t1045    Bacteria Bacteroide… Bact… Bact… Riken… Alis… putred… TACGGAG…
##  8 t1049 t1049    Bacteria Bacteroide… <NA>  <NA>  <NA>   <NA>  <NA>    TACGGAG…
##  9 t1089 t1089    Bacteria Proteobact… Delt… Desu… Desul… Bilo… wadswo… TACGGAG…
## 10 t1090 t1090    Bacteria Proteobact… Delt… Desu… Desul… Bilo… <NA>    TACGGAG…
## # ℹ 678 more rows
## # ℹ 2 more variables: genus_oldtaxonomy <chr>, genus_nosubgroups <chr>
## 
## attr(,"class")
## [1] "tidytacos"

Podemos visualizar la diversidad alpha de muestras de participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no:

vdata_rar %>%
  samples() %>%
  ggplot(aes(x = Sexual.Intercourse.24hours, y = observed, fill = Sexual.Intercourse.24hours)) +
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(height = NULL)

Análisis de coordenadas principales (PCoA)

Nos gustaría abordar las diferencias entre muestras de participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no. También estamos más interesados en los géneros que en los ASV.

Una PCoA podría ofrecer información:

vdata_genus <- vdata %>%
  aggregate_taxa(rank = "genus")

tacoplot_ord_ly(vdata_rar, Sexual.Intercourse.24hours, samplenames = sample, 
                dim = 2)
## Warning in as.dist.default(xdis): non-square matrix

La función aggregate_taxa fusiona todas las filas de la tabla de taxones en un nivel taxonómico específico, en este caso el nivel de género. Como ocurre con todas las funciones de tidytacos, todas las demás tablas del objeto tidytacos se ajustan en consecuencia.

La función tacoplot_ord_ly determinará la abundancia relativa de taxones en las muestras y luego utilizará las diferencias de Bray-Curtis para ordenar muestras en un espacio bidimensional (o tridimensional) según su composición taxonómica. La adición argumental de plotly “_ly” hace que la figura sea interactiva, lo cual es realmente bueno para el trabajo exploratorio. Esto también funciona para otras funciones de visualización.

Relación entre la composición de la comunidad microbiana y variables.

La siguiente pregunta lógica es hasta qué punto la actividad sexual reciente (dentro de las 24 horas) determina la variabilidad de la composición de la comunidad microbiana. No olvidemos que la composición microbiana vaginal varia con la edad, por eso incluyamos edad en el modelo.

perform_adonis(vdata_genus, c("General.Age", "Sexual.Intercourse.24hours"))
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = as.formula(paste("counts_matrix", formula_RHS, sep = " ~ ")), data = metadata, permutations = permutations)
##                             Df SumOfSqs      R2      F Pr(>F)    
## General.Age                  1    1.717 0.02676 5.5016  0.001 ***
## Sexual.Intercourse.24hours   1    0.956 0.01489 3.0607  0.015 *  
## Residual                   197   61.500 0.95835                  
## Total                      199   64.173 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La función perform_adonis realizará un ANOVA PERMutacional para determinar el efecto de las variables de ¨samples¨ en las diferencias de Bray-Curtis de las comunidades. El resultado muestra que la edad es un contribuyente a la composición de la comunidad microbiana (R cuadrado = 0.02676).

Análisis de abundancia diferencial

A continuación, nos gustaría saber cuáles de los 20 géneros más abundantes son significativamente más abundantes en las participantes que realizaron actividad sexual dentro de las 24 horas antes de tomada las muestras versus las que no.

vdata_genus <- vdata_genus %>% add_codifab(Sexual.Intercourse.24hours, 
                                           max_taxa = 20)
vdata_genus$taxon_pairs <- filter(vdata_genus$taxon_pairs, wilcox_p < 0.05)
tacoplot_codifab(vdata_genus, FALSE_vs_TRUE)

La función add_codifab agregará una tabla llamada taxon_pairs al objeto tidytacos, con la abundancia diferencial del taxón entre las dos condiciones (con respecto al taxón de referencia), para cada par de un taxón y un taxón de referencia.

La función tacoplot_codifab devuelve un gráfico para visualizar la abundancia diferencial de taxones entre condiciones, en comparación con todos los demás taxones como referencia. Podemos observar que es más probable que Staphiloccus y L. iners group sea típico encontrar cuando una participante ha tenido actividad sexual dentro de las 24 horas antes de tomada las muestras.

Es de destacar que existen muchos métodos de análisis de abundancia diferencial y ninguno de ellos es perfecto. Interprete tus resultados con cuidado.


Descargo de responsabilidad: Esta guía es una adaptación de la versión en inglés del tutorial de tidytacos.